function [K, grads] = get_K(x,xb, p,pb, theta, k0, c)
    % x \in R^{4xN}, spacetime convention: (ct, x, y, z).四维时空
    % p \in R^{3xN}
    % theta = struct array{rho, v, mu}: {model weight, speed, concentration} 

    [~, N] = size(x);
    [~, Nb] = size(xb);
    K = zeros(Nb, N); 
    grads = struct('rho', {}, 'v', {}, 'mu', {});

    for idx_row = 1:Nb
        for idx_col = 1:N

            s1 = xb(2:4, idx_row);   s2 = x(2:4, 33);   % Extract spacetime coordinates. 
            t1 = xb(1, idx_row)/c;   t2 = x(1, 33)/c;
            p1 = pb(:, idx_row);     p2 = p(:, idx_col);
            
            rho = theta.rho;        % Channel power, positive scalar
            v   = theta.v;          % Velocity      R^3, the derivative is also R^3. 
            mu  = theta.mu;         % Concentration R^3, the derivative is also R^3.
            mu_len = norm(mu);

            z   = (s1-s2) + v*(t1-t2) - 1i*mu/k0;   % Complex antenna displacement vector. 
% z   = -10*6.25e-4- 1i*mu/k0;
% z   = - 1i*mu/k0;
            [R, dR] = get_R(k0*z);                  % Compute 3x3 correlation matrix R(x1, x2). 
            [C, dC] = get_C(mu_len);              
            pRp = p1.'*R*p2;                        % Compute the projection of R onto polarizations p1, p2. 
            
            rho*(pRp)/C

            K(idx_row, idx_col) = rho*(pRp)/C;%？
            % K(idx_row, idx_col) = rho*pRp;          

            % dK_rho = pRp/sinc(1i*mu_len/pi);
            dK_rho = pRp; 
            dK_mu = zeros(3, 1);
            dK_v  = zeros(3, 1);
            for idx = 1:3
                dK_mu(idx) = rho*(-1i/C*(p1.'*dR{idx}*p2)-(dC*mu(idx))/(C^2*mu_len)*pRp);
                % dK_mu(idx) = rho*(-1i*(p1.'*dR{idx}*p2));   % Compute the derivative of K w.r.t. mu(1,2,3). 
                dK_v(idx) = rho/C*k0*(t2-t1)*(p1.'*dR{idx}*p2);
                % dK_v(idx) = rho*k0*(t2-t1)*(p1.'*dR{idx}*p2); 
            end

            grads(idx_row, idx_col) = struct('rho', dK_rho, 'v', dK_v, 'mu', dK_mu);
        end
    end

end



